Detection of periodic gravitational wave sources by 
Hough transform in the / vs / plane 



r.Antonucci, P.Astone, S.D'Antonio, S.Prasca, C.Palomba 
INFN Rome "La Sapienza", INFN Rome 2 "Tor Vergata" 
and University "La Sapienza", Rome 



Abstract. In the hierarchical search for periodic sources of gravitational waves, 
the candidate selection, in the incoherent step, can be performed with Hough 
transform procedures. In this paper we analyze the problem of sensitivity loss 
due to discretization of the parameters space vs computing cost, comparing the 
properties of the sky Hough procedure with those of a new frequency Hough, 
which is based on a transformation from the time - observed frequency plane 
to the source frequency - spin down plane. Results on simulated peakmaps 
suggest various advantages in favor of the use of the frequency Hough. The ones 
which show up to really make the difference are 1) the possibility to enhance the 
frequency resolution without relevantly affecting the computing cost. This reduces 
the digitization effects; 2) the excess of candidates due to local disturbances in 
some places of the sky map. They do not affect the new analysis because each 
map is constructed for only one position in the sky. 
Pacs. numbers: 04.80Nn,07.05Kf,97.60Jd 



1. Introduction 

One of the main goals of the search for periodic isolated sources of gravitational 
waves (g.w.) is to perform all sky surveys, based on "blind searches", where the 
source parameters are unknown. In this case hierarchical procedures are applied, 
based on a sequence of increasing resolution steps. In this paper we study in details 
the problem of sensitivity loss due to discretization of parameters and to the needs 
to limit the computing cost, with Hough procedures. In particular, we propose and 
study the characteristics of a frequency Hough procedure, designed mainly to reduce 
the discretization problem, and we compare it with the sky Hough procedure, which 
is actually used in the Virgo collaboration. 

The paper is organized as follows: in Sect. 2 we present the basic scheme of the Rome 
hierarchical procedure, based on the main idea of coincidences among subsets of data; 
in Sect. 3 we discuss the limits due to digitization of the sky Hough procedure; 
in Sects. 4, 5 we present the new frequency Hough procedure, discussing details 
its implementation and its basic characteristics; in Sect. 6 we present the study of 
amplitude losses due to digitization, and thus efficiencies, for both the procedures. 
Conclusions and comments are given in Sect. 7. 
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2. Scheme of the hierarchical procedure and Hough transforms 

Hierarchical procedures, based on Hough transform algorithms, are applied by various 
groups in the g.w. community. See, for example, references [HE]. There are various 
ways of implementing the hierarchical procedure and the Hough transform. 

2.1. Hough transform basic characteristics 

The Hough transform is a linear transform that is used to recognize the parameters of 
the analytical description of a curve from the position of some points on it. It operates 
on an "image" of points, in our case the peakmap in the time-frequency plane. For 
each peak of this map we increase a set of bins of a multi-dimensional histogram (in 
our case a two-dimensional histogram) defined on the parameters space, called the 
Hough map. In the old procedure, the parameter space was the position of the source, 
i.e. the celestial sphere, and we fixed the spin down value for each Hough map. In 
the new one, the parameter space is the plane / — /, and for each Hough map, we fix 
the position of the source. The mapping (i.e. which points of the Hough map must be 
increased for a certain point in the peakmap) can be done in different ways: we use 
always what we call the "biuni vocal mapping" , i.e. a mapping in which every point 
in the Hough map derive from a single point of the peakmap at a given time. It is 
easy to demonstrate that in this case the mapping is also uniform, i.e. in the case of 
uniformly distributed random dots in the peakmap, the expected value of the Hough 
map is a constant (for all parameter value). This value, depending on the number 
N of the spectra of the peakmap and on the mapping, defines the "noise" of the map. 
It is binomially distributed with parameters N and p = ^./N . 

2. 2. Scheme of the hierarchical procedure 

We will refer here to the Rome scheme, presently used in Virgo data. Fig. [1] shows 
the basic scheme of the Rome hierarchical procedure. Details on the main aspects of 
the procedure are given in references [31 l4l [5] . After data cleaning (short time domain 
disturbances removal) and "Short FFTs Data Base" (SFDB) creation, peakmaps are 
computed, using a very refined auto-regressive algorithm to equalize the spectral data 
by an appropriate follow-up of the noise. Peakmaps are frequency vs time maps, 
obtained from equalized spectra by selecting all the local maxima above a chosen 
threshold. An accurate cleaning of peakmaps, by removing known noise lines and the 
more persistent lines, is needed and its implementation is critical for the next step 
analysis. On the cleaned peakmaps, methods of peaks detection are applied. That is, 
transformation from the input plane to the Hough plane, thresholding and first order 
candidates selection. Candidate parameters are defined by source frequency, celestial 
coordinates, first spin-down parameter. The need for coincidences among candidates 
obtained in different subsets of data (two in the scheme of Fig. [1]) has been discussed 
in references O [7l [8] . This method is very efficient to reduce the number of spurious 
candidates at a fixed threshold. Thus, for a given false alarm probability, we can 
lower the threshold -with respect to the choice of not doing coincidences- gaining in 
detection efficiency. The method has a better efficiency when the data sets have similar 
sensitivities. After the coincidence, the survived candidates are analyzed coherently 
with longer FFTs on corrected data. Then the spectral filtering is used to take into 
account the spread of the power in five bands, as explained in reference [5]. Finally, 
second order candidates are produced. 
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Figure 1. The basic scheme of the hierarchical procedure, used in the Virgo 
collaboration 

3. The problem 

As stated before, the sky Hough method shows amphtude losses, and thus loss of 
sensitivity, which are due to digitization of parameters. This effect shows up mainly 
for the complexity of the transform together with the need of reducing the computing 
cost: 

• the method is based on a transform between the time-frequency peakmap and 
the celestial sphere. It is not simple for the non linearity of the mapping; 

• to reduce the computational effort, we need to use "look-up tables" which 
introduce further digitization errors; 

• to reduce the computational effort, fast algorithms have been developed, which 
require the use of a rectangular grid to map the sky. Compared to the "optimal" 
(see later) grid, the rectangular one has over-resolution in some regions of the 
sky. This leads also to a higher number of candidates. 

• the use of the celestial map as the space to spot the candidates is very prone to 
artifacts, see [B]: some regions are always "privileged", that is they have a higher 
candidates number with respect to the expectation. The problem arises because 
each Hough map is constructed over the whole sky. 

Hence, it seemed important the study of alternative procedures. Given the observation 
that most of the problems are related to the complexity of the transformation, we 
exploit the possibility of the use of a different but simpler transformation. A part 
the simplicity of the transformation we obviously need to study a procedure which 
is less, or equivalently, computationally expensive. Therefore we studied a procedure 
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which has a better, or equivalent, sensitivity, at the same computational cost of the 
sky Hough. 



4. The frequency Hough 

The transformation we propose transforms the Time - observed frequency plane 
into the Source frequency - Spin down plane. Let' s go into details. If / is 

the frequency (Doppler corrected for a given sky direction), /o the source intrinsic 
frequency, d = f the first spin-down parameter, t the time at the detector and to a 
reference time, we have that 

f = fo + d{t-to) (1) 
a straight line in the Hough plane. We then get the following: 

d= + (2) 

t — to t — to 

Each point in the input plane {t — to,f), that is a peak in the Doppler shifted peakmap, 
is transformed into a straight line in the Hough (/o, d) plane, with slope —l/{t — to)- 
The slope depends on the choice of the reference time. If we choose to equal to the 
beginning time of the data we analyze, then the slope is always negative and inversely 
proportional to the time gap. 

This is the choice we have done here. In addition, considering the width A/ of the 
frequency bins in the input plane we notice that each peak is transformed into a stripe 
among two parallel straight lines 

t — to t — to t — to t — to 

It is a linear transformation. Now the input plane is obtained from the original 
peakmap by correcting it for the Doppler shift due to the Earth revolution and 
rotation, for each point in the sky grid we need to analyze. Thus "time" is the time 
at the detector and "frequency" the observed frequency, after the Doppler correction. 
But, as each SFDB is short enough to not be affected by a time- varying Doppler 
shift, then the Doppler effect removal from the original peakmap, obtained from the 
collection of all the SFDB data, reduces to a very simple "shifting" procedure of the 
peakmap bins. In the analysis scheme, this bins shift is part of the Hough procedure. 
In the following, we give details on the construction of the map. 



4-1- The direct differential method 

The frequency Hough map is constructed using the "direct differential method" , as is 
done with the sky Hough. With this method, instead of building directly the Hough 
map, one builds a map that, if "integrated" (i.e. summed over bins from left to right), 
gives the Hough map. This is important to minimize the number of floating point 
operations. As already explained, for each sky position, the input peakmap is got 
from the original one by shifting bins to correct for the Doppler effect. The sky is 
sampled with a non uniform covering grid, which will be later discussed. Here we 
explain in detail the technique, by giving the sequence of operations: 

• for each point in the sky grid, for each coordinate in the input plane {t — to, f) 
and for each spin- down value d, 
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the map is incremented by 1 in the point 



fo = f 




to) 



and decremented by 1 in the point 



/o = / 



+ 




to) 



Hence, for each sky position, a differential map is constructed. The sum of 
the bins along the frequency direction is then performed to construct the final 
integral map. This two dimensional histogram is the frequency Hough map. In the 
algorithm implementation we plan to divide the input peakmap into 10 Hz bands, 
thus constructing, for each position in the sky, a different Hough map every 10 Hz. 
In case there is the need to exploit higher order one spin down parameters, one (or 
more) loop(s) has (have) to be added to the sequence of operations, to scan the discrete 
set of values of the new paramctcr(s). This clearly influences the computing cost, but 
does not change the basics of the method. 

5. Main characteristics: frequency resolution and sky grid 

Let's first discuss two peculiar aspects of this new method, which are the basis of its 
appeal. 

5.1. Increasing the frequency resolution 

From the given analysis scheme, it is easy to see that the frequency resolution for the 
estimation of the somxic frequency /o can be enhanced, with respect to the binning 
frequency A/, without relevantly affecting the computational effort. In fact, the use 
of a resolution 



with r >= 1, affects only the size of the Hough map. This has a computational cost 

only when summing over the bins to construct the integral map from the differential 
one. But we notice that the total cost of the construction of the Hough map is due 
to the construction of the differential map, dominated by the number of peaks in the 
peakmap and to the construction of the integral map, dominated by the number of 
bins. The former, in all practical cases, is the one which dominates. 
The possibility to enhance the frequency resolution results to be, as will be shown 
in the next sections, a very important peculiarity of the new method. It which 
enhances considerably the efficiency, by reducing the digitalization effect. The same 
in the sky Hough procedure would have a relevant computational cost. Regarding 
the increasing of the spin down resolution, it would cost for both the procedures: the 
better the resolution in the spin down estimation the higher is the number of loops of 
the procedures. 

5.2. The grid on the sky 

Here we describe how we construct the grid on the sky. Suppose two sources, at the 
same frequency fo and same latitude /3. Their angular delay 7 with respect to the 



A/o = Af/r 



(4) 
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detector rotation produces a time delay At — j /Qaj-i,. The two sources will then have 
the same frequency variation at the detector, which is the classical equation due to 
the Doppler effect, 

. J, „ Vorb COS (3 

— Jo 

c 

but with time delay At. The observed frequency difference has thus a maximum value 
which is given by 

A/ = 2/0 (5) 

c 

Thus the angular resolution is, in radians: 

7 = l/(7Vz5 cos/3) (6) 

where Nd is the number of points in the Doppler band for a signal of max frequency 
/o: 

ND = 2foVorb/{cAf) (7) 

and Vorb I c ~ lO""*. 

We now repeat the same reasoning, supposing the two sources, at the same frequency 
/o and same longitude A. The two sources will have the same frequency variation at 
the detector, now given by but with an angular delay 7'. The observed frequency 
difference has a maximum value which is: 

^^^^^W^ (8) 
c 

We obtain for the angular resolution, in radians: 

7' = 1/{Nd sin/3) (9) 
Using eqs. |6] and [9] we get: 

AA = 1/(A^D cos/3) (10) 

A/3 = l/(iVi3 sin/3) (11) 

Using these equations we construct the grid on the sky, which we call the "optimal" 
grid. The points of the grid are not uniformly distributed. With a simulation, we have 
estimated the the number of points in the grid Ngi-y, which is, in the high frequency 
limit: 

N sky— K sky (12) 

TT 

Ksky is an extra resolution factor, which can be greater than 1, to enhance the 
efficiency, but even less than 1, to save computing cost, obviously worsening the 
efficiency. FigH] shows the optimal sky grid, for Njj = 20 (which corresponds to a 
source frequency /o = 100 Hz). As already said, the grid used in the sky Hough 
method, is not optimal, but rectangular, to use fastest computing algorithms. The 
number of points in this rectangular grid is: 

Nsky=Ksky27r^Nl (13) 

which is, asymptotically, a factor tt higher then the number of points of the optimal 
grid. In fact this grid has to be over resolved to maintain the same sensitivity of the 
corresponding optimal grid. Further, we note that this over resolution produces a 
higher number of candidates from certain sky positions. 
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Figure 2. Optimal sky map, for Nd ~ 20, Ksky = 1- x-axis: Ecliptical 
longitude, degrees, from to 400; y-axis: Ecliptical latitude, degrees, from 
-100 to 100; The number of points in the map is Nsky=2902. 



5. 3. Effect of artifacts 

The sensitivity of the sky Hough procedure is affected by artifacts, i.e. an excess of 
candidates in some places of the sky map, which are due to local spectral disturbances. 
The effect can't be eliminated because each map is constructed over the whole sky, 
and hence the threshold for candidate selection has to be the same for the whole 
sky. Using the frequency Hough procedure this effect disappears because each map 
is constructed for only one position in the sky. So, because of the adaptivity of the 
threshold, if a sky region gives an excess of candidates, the threshold is raised and 
then there is a loss in sensitivity only for that sky region. 

6. Study the efficiency of the methods 

We are now ready to enter into details by studying the efhciency of both the methods, 
by the use of simulations. Figure [3] is an example of how a frequency Hough map 
looks like, having injected into white noise three signals, at different frequencies and 
spin-down. 

6.1. Loss in amplitude vs frequency over resolution factor 

To study the efficiency of the methods, as a function of the frequency over resolution 
factor, we have simulated a signal in the absence of noise. The reason for this is 
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f.. [Hz] 

Figure 3. An example of a frequency Hough map, having injected into white 
noise three signals, at different frequencies and spin-down, x-axis: Estimated 
source frequency, Hz, from to 100; y-axis: Spin down, Hz/s, from 0.02 to 
-0.02. 



that we were interested in studying only the losses due to the discretization errors. 
The parameters chosen for the simulation are similar to actual situations (detector 
parameters, source expected parameters). The parameters of the simulation are shown 
in Table[T] Fig. |3]shows the amplitude loss versus the frequency over resolution factor 
r. The loss was calculated as the average value of all the peaks found in the 500 
spectra (it is important to remember that our procedure considers peaks only the 
maxima above threshold). The result is clear: using r = 10 the amplitude loss is 
3.6 % (the efficiency 96.4%), while with r = 1, which is the only practically possible 
choice of the sky Hough, the amplitude loss is 11.6 % (the efficiency 88.4%). From the 
figure, we notice that there is no further gain of increasing the over resolution factor 
over 10. Thus, we fixed to 10 the over resolution factor for the frequency Hough. In 
next simulations, results with r = 10 are thus for the frequency Hough, results with 
r = 1 are for the sky Hough. 

6.2. Loss in amplitude vs the spin down resolution 

Once we have fixed the frequency over resolution factor we wanted to quantify how 
the increasing of the spin down resolution from the nominal one would affect the 
sensitivity. The results are in Fig. [5l which shows the loss in amplitude vs the spin 
down over resolution factor, for both the cases r = I, sky Hough, and r = 10, frequency 
Hough. It can be noticed that, in the case of the frequency Hough, even for the worst 
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Figure 4. Loss in amplitude vs the frequency over-resolution factor. 
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Figure 5. Loss in amplitude vs the spin down over-resolution factor, for the 
two cases: r = 10 (left) and r = 1 (right). 



analyzed situation, which corresponds to the nominal spin down step 

the loss is quite small. Is is 3.6 % (the efSciency 96.4%). The situation is worst for 
the sky Hough, where the loss in amplitude at the nominal spin down step is 11.6 % 
(the efficiency 88.4%). The improvement obtained by a better spin down resolution 
is not so important, as can be seen from the figure. It seems reasonable, given the 
observation that increasing the spin down resolution has a computational cost for both 
the methods, to use the nominal Ad resolution (x-axis equal to 1 in the figure). 



6.3. Loss in amplitude vs sky grid resolution 

To study the loss due to the sky grid resolution, we have simulated 50 signals, ran- 
domly distributed over the sky. We have then looked for results using the optimal 
grid, again registering the average value of all the detected peaks. In what follows, we 
suppose to use the optimal grid for both the procedures, sky and frequency Hough. 
Fig[6] shows the amplitude losses, as a function of the over resolution sky map factor 
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Figure 6. Loss in amplitude versus the sl<y map resolution factor Ksky The 
figures compare the loss when r = 10 (left), frequency Hough, and when r = 1 
(right), sky Hough. 



A/ [Hz] 


Source frequency [Hz] 


d [Hz/s] 


Ad [Hz/s] 


10-3 


50 


5-10-8 


2 - 10-9 



Table 1. Parameters used in the simulation. The number of SFDBs is 500. The 
nominal spin down resolution is Ad = . 





Sky 


Frequency and spin-down 


Total 


Sky Hough method 


0.860 


0.88 


0.757 (squared 0.573) 


Frequency Hough method 


0.900 


0.964 


0.868 (squared 0.753) 



Table 2. Comparison of h efficiencies. The table reports, for both the methods, 
the partial (for unknown sky position, and for unknown frequency and spin down 
values), and the total efficiencies. 



Ksky, in the two cases of r = 10 (left), frequency Hough, and r = 1 (right), sky Hough. 
The amplitude loss, for Ksky — 1, is 10% for the frequency Hough, and 14%, for the 
sky Hough. Again, a better efficiency for the new procedure. We notice that the use 
of an over resolution for the sky map, would have an impact on the computing cost, 
with both the procedures. 



7. Summary of the results and conclusions 

Results of this study are summarized in Table [2] 
We see that the ratio of the amplitude efficiencies is 

freq. Hough ^ ^ 
sky Hough 
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which in power is 1.317. Prom this, we can compute the gain in computing cost for the 
same sensitivity. Let us firstly recall that the h sensitivity in the hierarchical search is 
proportional to Tcoh* , and the computing cost to T^oh ■ Thus, the "equivalent FFT" 
length factor is 1.148^=1.734 and the gain in computing cost is 1.734^=5.2 (that is, 
the ratio of computing costs needed to have the same h sensitivity). 

Let's conclude with two further considerations, regarding the characteristics of 
the frequency Hough procedure: 

• the adaptivity, that is the weight of peaks to consider the noise level and the gain 

due to the antenna pattern toward a direction, is, with this approach, immediate 
and very simple, as each Hough map is done for a single sky position, it has 
been shown, with the sky Hough, that the adaptivity of the procedure is a very 
important task for the analysis; 

• this new procedure is appropriate also for all those situations in which the source 
position is known and we should estimate only source frequency and spin down; 

• with a proper choice of parameters, it is also possible to detect and hence remove 
spurious signals, with a constant or linearly varying frequency. 

On the latter point, wc arc now working to study the efficiency of this method in terms 
of rejection of spurious lines in the peakmap. We know that this is a very critical task 
for the analysis, since the presence of spurious lines highly affects the sensitivity of the 
search. Wc expect this new method to be much more insensitive to the presence of 
spurious lines, since in the chosen Hough plane spurious lines and g.w. signals should 
have a very different and well separable behavior. 
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